Dynamic Supercoiling Bifurcations of 
Growing Elastic Filaments 



Charles W. Wolgemuth ^ 

^Department of Cell Biology, University of Connecticut Health Center, 

Farmington, CT 06030 

Raymond E. Goldstein^ 

^Department of Physics and Program in Applied Mathematics, University of 

Arizona, Tucson, AZ 85721 

Thomas R. Powers ^ 

^Division of Engineering, Brown University, Providence, RI 02912 



Abstract 

Certain bacteria form filamentous colonies when the cells fail to separate after di- 
viding, bi Bacillus subtilis, Bacillus thermus, and cyanobacteria, the filaments can 
wrap into complex supercoiled structures as the cells grow. The structures may be 
solenoids or plectonemes, with or without branches in the latter case. Any micro- 
scopic theory of these morphological instabilities must address the nature of pattern 
selection in the presence of growth, for growth renders the problem nonautonomous 
and the bifurcations dynamic. To gain insight into these phenomena, we formulate 
a general theory for growing elastic filaments with bending and twisting resistance 
in a viscous medium, and study an illustrative model problem: a growing filament 
with preferred twist, closed into a loop. Growth depletes the twist, inducing a twist 
strain. The closure of the loop prevents the filament from unwinding back to the 
preferred twist; instead, twist relaxation is accomplished by the formation of su- 
percoils. Growth also produces viscous stresses on the filament which even in the 
absence of twist produce buckling instabilities. Our linear stability analysis and nu- 
merical studies reveal two dynamic regimes. For small intrinsic twist the instability 
is akin to Euler buckling, leading to solenoidal structures, while for large twist it is 
like the classic writhing of a twisted filament, producing plectonemic windings. This 
model may apply to situations in which supercoils form only, or more readily, when 
axial rotation of filaments is blocked. Applications to specific biological systems are 
proposed. 
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1 Introduction 



Many structures in the biological world grow so slowly that they adopt a 
shape that can be considered as a minimizer of some configurational energy 
associated solely with the internal structure. The logarithmic spiral of the 
nautilus shell is an example. It enlarges through a process of differential growth 
whereby its shape represents the accumulated history of identical events, save 
for scale changes [1]. Thus, the microscopic rules of growth are essentially 
imchanging as the three-dimensional form develops, and the properties of the 
external environment do not fundamentally determine the form. 

Environmental effects on the development of biological forms are well known. 

Consider the formation of tendril perversions in climbing vines [2,3]. A per- 
version is a junction between regions of opposite helix handedness that forms 
as an initially straight tendril first attaches to a support structure and then, 
through the activation of tension-sensitive receptors, undergoes a helical in- 
stabihty. The constraint of fixed ends enforces the formation of a structure 
with zero net twist, consisting of concatenated regions of opposite chirality 
joined by a transition region — the perversion. Thus, the interaction with the 
external environment fundamentally alters the pattern formation through the 
tension induced by a point contact. 

Even more complex phenomena occur for those structures whose very process 
of growth induces an external force, as with motion through an environment 
that offers viscous resistance. These new forces can dramatically alter the 
ultimate configuration, and the present paper is a case study in such phenom- 
ena. We focus on a system in which the formation of patterns occurs through 
a finite-wavelength instability in which the process of growth introduces an 
intrinsic time dependence to the control parameters. The resulting bifurca- 
tion problem is nonautonomous and can exhibit a rich phenomenology as the 
intrinsic time scale of growth competes with those of the various modes of 
instability. This competition places the problem among the class of so-called 
"dynamic bifurcations," of which many examples are of continuing interest. 
These include instabilities in directional solidification in which the initial ac- 
celeration of the interface from rest provides the nonautonomous character [4] , 
fingering instabilities of magnetic fiuids under the infiuence of time-dependent 
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Fig. 1. Phase contrast micrographs of growing B. subtilis spores with attachment of 
both cell poles to the spore coat. The series consists of different filaments at various 
stages of growth. Scale bar = 10 fim. Figure courtesy of N. Mendelson; see [18]. 

magnetic fields [5,6], and separatrix crossing in Hamiltonian systems viewed 
as models for stellar fission [7]. 

Motivation for this focus comes from the phenomenon of supercoiling exhib- 
ited by filamentous colonies of Bacillus subtilis and other bacteria. B. subtilis 
cells are rod-shaped bacteria, typically four microns in length and slightly 
over one-half micron in diameter. Wild-type rod-shaped bacteria grow by ex- 
tending along the cylindrical axis of symmetry, and then dividing and sep- 
arating in the middle [8,9]. Under certain conditions, the cells of some mu- 
tant forms fail to separate upon replication, leading to a long chain of cells. 
Other species have also been observed to form chains, including Escherichia 
coli [10], Cyanobacteria [11], Myxococcus xanthus [12], and Mycobacterium tu- 
berculosis [13]. Under certain growth conditions, strains of B. subtilis, Bacil- 
lus stearothermophilus [14], Thermus [15], and Mastigocladus laminosus [16] 
form complex braided structures. Helical (or, in the jargon of DNA bio- 
physics, "solenoidal" ) morphologies have also been observed [17]. Of these 
examples, the supercoiled structures of B. subtilis have been studied most 
extensively [18,19]. Recent experimental work [20,21] has indicated that the 
cell wall of B. subtilis contains helical protein structures. These may supply 
the molecular imprinting responsible for this morphological development, in 
a manner analogous to the way microtubules control macroscopic handedness 
in certain plants [22], but as of yet no successful microscopic theory for the 
formation of these supercoiled structures exists. 

Throughout much of the development of complex structures in B. subtilis, the 
length L of the elongating chain, a single cell thick, grows exponentially in 
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time, L oc exp((Tt), with the growth rate a ~ 2 x 10~^ s~^. Moreover, material 
cross-sections of the cells rotate relative to each other [23], so that the angle 
describing the relative orientation of any two material cross-sections is pro- 
portional to their exponentially increasing separation. In one of the earliest 
observations of supercoiling in B. subtilis, the ends of the fiber adhered to a 
spore coat, prohibiting axial rotation of the ends [18]. Mendelson supposed 
that this blocked rotation leads to torsional stress, eventually causing super- 
coiling [18] (See Fig. 1). Since this discovery, it has been shown that adhesion 
is not required for supercoiling [19], and that many different factors such as 
temperature, pH, and the concentration of ions such as magnesium and am- 
monium affect the morphology of the coils [24,25]. The evolution of the coils 
after the formation of the first braid is remarkable. The first plectonemic braid, 
essentially a filament that is two cells thick, continues to grow and eventually 
reaches a critical length of order 100 /xm, after which it supercoils to form 
another braided structure which is four cells thick. This process can repeat 
many times, leading to a hierarchy of braids and eventually a macroscopic 
object. 

In this paper, we study in detail Mendelson's blocked rotation mechanism of 
supercoiling in a growing closed loop, and focus on the formation of the initial 
braid. While not yet explaining the microscopic physical origin of the coiling 
instability, we elucidate the rich dynamics that occurs when growth competes 
with blocked rotation, and thereby help constrain more detailed theories. 

With this focus we exploit several simplifying assumptions in our analysis. 
First, the chain of cells is treated as an elastic filament with uniform properties 
along its length. There is evidence that this assumption holds until times 
comparable to a few doubling times (~ o"^^), but may be violated later. For 
example, when the chirality of nutrient molecules in the growth medium is 
reversed, the supercoils unwind and even begin to wrap up in the opposite 
handedness, but the hairpin bends from the original braid remain [26] . Thus, 
some of the deformation of the growing filament becomes permanent. This 
phenomenon is reminiscent of the morphological development of plant tendrils, 
in which young and fiexible tendrils age with time, becoming woody and locked 
in a fixed shape [27]. On the time scales that will concern us, t << a~^, single 
fibers have been shown to behave like elastic rods, with a bending modulus 
A — 10~^^ erg-cm [28], as discussed further in section 2.4. 

A second major simplification we introduce is to treat the growth rate as 
constant in time, independent of stress and filament geometry. Again, the 
permanent hairpin bends in the chirality reversal experiment show that this 
assumption cannot hold everywhere along the filament for all times after the 
first braid forms. 

Together with the observation that viscous effects dominate inertial effects in 
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the low Reynolds number environment of the growing fibers, these assump- 
tions lead to the model studied in this paper: an elastic ring with intrinsic 
twist suspended in a viscous fluid and lengthening at an exponential rate. Our 
work is complementary to that of other investigators. For example, Shelley and 
Ueda studied the Euler-like buckling of a growing liquid crystal filament, using 
a local drag model for the linear stability of a growing loop [29], and incorpo- 
rating nonlocal Stokesian hydrodynamics to study the pattern formation [30]. 
Drasdo studied similar patterns in the context of the growth of single-cell-layer 
tissue sheets [31]. Klapper has studied inertial writhing instabilities of open 
rods subject to exponential growth, as well as the relaxation to equilibrium 
of twisted rings in the absence of growth [32]. Goriely and Tabor [33] intro- 
duced the idea of twist depletion as a possible mechanism driving buckling in 
B. suhtilis. 

Our analysis begins with a generalization of the kinematics and dynamics of 
slender filaments to account for growth. Section 3 treats the growing elastic 
loop, beginning with a qualitative discussion of the instability. The linear 
stability analysis is greatly simplified by the use of the natural frame, so 
in this section we include a self-contained summary of the properties of the 
natural frame. We present a quasi- analytic treatment of the linear stability of 
the loop, and then present numerical simulations of the full nonlinear problem. 
Section 4 is the conclusion. 



2 Kinematics and Dynamics of Growing Rods 

2.1 Centerline Kinematics 

In this section we extend the standard kinematics and dynamics of elastic rods 
to allow for growth. Let s denote the arclength measured at time t from one 
end of an open rod, or from a fixed material point for a rod closed to form 
a loop. Then r(s, t) is the position in space of the centerline of the rod with 
arclength coordinate s at time t. Since material points on the rod centerline 
are convected along the rod by growth, fixed values of s do not correspond to 
fixed material points. We choose to label the material points of the centerline 
at all times by the arclength parameterization sq at a fixed time t = 0. We 
will study exponential growth, for which s — exp{at)so- Note that the partial 
derivatives with respect to s and t commute because s and t are independent 
variables: 



however, d/dt\sQ and d/ds\t do not commute since 



d_ 

di 



so 



d_ 
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ds 



so 



d_ 

ds 



(2) 



by the chain rule. 

The velocity of a material point is the time derivative of position at fixed Sq, 



dr 



(3) 



so 



It is convenient to develop the equations of motion in terms of s rather than 
So, so that 



dr 



ds 



so 



dr 

ds 



(4) 



To avoid confusion, we will explicitly denote which variables are fixed when 
finding the partial derivatives with respect to t. However, since partial deriva- 
tives with respect to s will always be taken at fixed t, we will write d/ds\t = 
d/ds. 

The first term of Eq. (4) corresponds to the velocity of the centerline in the 
absence of growth, and the second term arises from growth-induced convection. 
Henceforth we will specialize to exponential growth, for which ds/dt\sQ = crs 
and 



dr 

dt 



dr dr 

ds dt 



+ saes, 



(5) 



where 63 = dr/ds is the unit tangent vector of the centerline. Although the 
velocity of a material point may seem from Eq. (5) to depend on the arbitrary 
choice of origin for s, Eq. (3) shows that the velocity is manifestly independent 
of this choice. Note that Eq. (5) imphes that es- dv/ds — a. 



2.2 Choice of Frame and Growth Model 



In Kirchhoff rod theory [34], the configuration of a rod is completely specified 
by the orientation of the material orthonormal frame {ei, 62, 63}. The vectors 
ei and 62 point to material points on the rod surface (Fig. 2). As the rod 
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Fig. 2. The material frame. The unit vector is tangent to the rod centerline. The 
other members of the orthonormal frame, ei and 62, point to material points on the 
rod's surface. 



bends and twists, the positions of these material points change, causing the 
material frames to rotate [35]: 



9ej 

ds 



(6) 



dej 
dt 



(7) 



The vector Q describes the bending and twisting strain at a given instant, 
and (jj is the angular velocity of a material frame at a given material point 
Sq. In general, fl and uj depend on the choice of material frame. Once the 
choice is made for a given configuration, say the stress-free state, then the 
choice is specified for all configurations. In the classical rod theory without 
growth, it is natural to align ei and 62 with the principal axes of the cross- 
section. If the cross-section is circular, as we henceforth assume, and if the 
rod is straight in the absence of stress, then there are many equivalent natural 
choices. For example, if the rod aligns along the 2;-axis when it is stress-free, 
then {61,62} = {x, y} is natural. Any uniform rotation of this frame about z 
is equally convenient; all these choices lead to f2 = in the absence of stress. 

If the rod is curved in the stress-free state, then the direction of curvature 
breaks the rotational symmetry of the circular cross-section and provides a 
natural choice for the directions of ei and 62. E.g., if the rod has a helical 
stress-free state, then we may take ei = fi and §2 = h, where n and b are the 
unit normal and binormal of the Frenet-Serret frame [36,37,38]: 



^63 




ds 
dn 

ds 



= Kn 



^63 + rb 



(9) 
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Fig. 3. Two different modes of growth. In (a) and (b), the rod lengthens with- 
out rotation. In (c) and (d), cross-sections at fixed material points rotate with an 
ever-increasing angular velocity. 



dh 

ds 



-rn. 



(10) 



where k is the curvature and r the torsion. Thus, fl = ^062 + tqC^ for this 
frame in the absence of stress, where Kq is the curvature and tq is the torsion 
of the helix. 

When Kq = 0, the helix of the last example degenerates to a straight rod with 
spontaneous twist tq. We will show below that if = and the cross-section 
is circular, then tq can be eliminated from the equations of motion for an 
inextensible (non-growing) rod, and therefore does not affect the rod shape 
and dynamics. However, tq has physical meaning for a growing rod, even if 
hq = and the cross-section is circular. Figure 3 illustrates two kinematic 
possibilities for growth. Each sub-figure shows the stress-free configuration of 
a growing rod at two different times. In all cases, the left end of the rod has a 
fixed position and orientation. Two growth schemes are shown; in the growth 
scheme of Fig. 3a and 3b, the material frames are carried to greater values of 
z by growth and have no angular velocity. A line of material points parallel 
to the 2;-axis at time t remains parallel to the 2;-axis at time t + At (Fig. 3a). 
However, the pitch of a helical line of material points increases as the filament 
grows (Fig. 3b). Since the angular velocity is zero, the number of helical turns 
is constant. 



Figures 3c and 3d illustrate a growth scheme in which the cross-sections rotate 
with an angular velocity that increases with arclength. In this line 
of material points parallel to the z-axis at time t wraps around the rod at 
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time t + At (Fig. 3c). As the rod grows, the material frames are carried to 
greater vahies of z but also rotate relative to the fixed frame at s = 0. In this 
paper, we will study the growth model of Fig. 3c and d since it describes the 
relative rotation of the cross-sections of the B. subtilis fibers [39,40]. Define 
^(so) = cos~^(ei(0) • ei(so)) as the angle between the orientation of the (zero- 
stress) material frames at sq and s = 0. We will suppose that 9{sq) increases 
linearly with sq and exponentially in time with rate a: 

0{so) = ToSoexp{at) = TqS. (11) 

Since 9{so,t) and s{so,t) increase in time with the same exponential rate, a 
helical material line on the rod surface with pitch 2ti/tq remains a helical line 
with the same pitch as time passes (Fig. 3d). Thus, the natural choice for the 
material frame in the stress-free state is 



ei = cos(sro)x + sin(sro)y (12) 
62 = - sin(sro)x + cos(sro)y (13) 
§3 = 2. (14) 

Note that f2 = roes and uj = asTQe^ for this frame; oj ^ since the material 
frames must rotate to maintain zero stress as the rod grows (in contrast to the 
growth model of Fig. (a)). The parameter tq not only characterizes the con- 
figuration of the natural material frame (as in the inextensible helix example) 
but also the mode of growth. 

2.3 Compatibility relations 



Geometry relates the strain vector ft and the angular velocity u). To see how, 
consider Fig. 4. The lower curve represents the centerhne of the rod at time 
t, and the upper curve represents the centerline of the rod at time t + dt. 
The labelled points on the upper curve have the same material coordinate as 
the corresponding points on the lower curve; thus, the arclength parameter 
for p3 is 6 exp(adt) ~ s -|- sadt, while the arclength parameter for is {s + 
ds) exp((Tdt) (s -|- ds)(l -|- adt). Let Ri be the rotation matrix carrying 
the frame at pi to ps, /?2 the rotation matrix carrying the frame at to p4, 
R3 the rotation matrix carrying the frame at pi to p2, and /?4 the rotation 
matrix carrying the frame at p2 to ^4. Furthermore, let J and K denote the 
infinitesimal rotation matrices associated with the rotation vectors ft and uj 
respectively {e.g. Jap — ea/s-yft-y, where Ca/s-y is the alternating symbol). From 
the definitions of the rotation matrices, R2R1 — RaRs, or 



/ + (1 + adt)dsJ{s + sadt, t + dt)] [/ + dtK{s, t) 
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l + dtK{s + ds,t) I + dsJ{s, t) 



(15) 



where / is the identity matrix. Expanding Eq. (15) to C(dsdt), we find 



dJ 

+ as— + aJ+[J, K] 

OS 

s 

+ aJ+[J, K], 

so 



dK 


dJ 


ds 


dt 


dK 


dJ 


ds 


dt 



(16) 
(17) 



where [, ] is the commutator. In terms of components in the material frame, 



dQi 




dui 


dt 


so 


ds 






duj2 


dt 


so 


ds 


dQs 




dujs 


dt 


so 


ds 



— aVls + niUJ2 — ^2'^i- 



(18) 
(19) 
(20) 



These compatibility relations show how strain changes in time due to non- 
uniform rotation rates (the first term in each of Eqs. (18)-(20)), growth (the 
second term in each of Eqs. (18)-(20)), and the geometric couphng between 
twisting and bending (the last two terms in each of Eqs. (18)-(20)). The 
compatibiUty relation Eq. (20) will be used below to determine the dynamics 
of the twist strain Q^- Note that Eq. (20) can be re- written in a form valid for 
arbitrary growth laws [41,42,43] 



dn. 



dt 



so 



doj-i ^ dr dv dr d'^r dw 
— - - ^3 — ■ \ X ■ — . 

ds ds ds ds ds'^ ds 



(21) 



2.4 Constitutive relations 



The Kirchhoff constitutive relations for a rod with intrinsic twist relate the 
moment on a cross-section of the rod to the strain [44] : 

^^^^ 

Note that dr/ds x d'^v/ds'^ = f2iei -|- (^262. Equation (22) impfics that in the 
stress-free state (defined by M = 0), ei and 62 rotate around §3 with rate 
To — e2 • dei/ds, and 63 is constant. 
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4 t+dt 



t 



Fig. 4. Geometrical origin of compatibility relations. 

The force and moment balances for a growing rod are the same as the balances 
for an inextensible rod: 



dF 

OS 



dM 

ds 



+ 63 X F + m, 



cxt 



0, 



(23) 
(24) 



where F{s,t) is the force the internal elastic stresses exert through the cross 
section at s on the portion of the rod with arclength less then s. The external 
force per unit length fext and moment per unit length nicxt are measured per 
unit arclength. In an alternative but equivalent formulation, the elastic force 
and moment per unit length arise from variational derivatives of the energy 



E 



ds 



^ft;2_^2(fi_fiQ)2_A(s,t) 



(25) 



where = |(9^r/(9s^p = Ql + is the square of the curvature, and A is 
the Lagrange multiplier associated with the constraint of prescribed length, 
L = Loexp{at) [41,42]. 

The external force per unit length fext consists of a viscous drag force per 
unit length and an artificial short-ranged repulsive force that prevents self- 
crossing: fext = fvisc + fseif- The external moment per unit length is purely 
viscous: mext = rHvisc- The repulsive force takes the form 



Iself 



/3(r(s) - r{s')) 



ds'. 



\s-s'\>S 



ris 



ris' 



(26) 
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where 5 is a short-distance cutoff, n = 14 simulates the repulsive part of a 
Leonard- Jones potential, and i3=180 is sufficient to keep the filament from self 
crossing. The viscous force per unit length fyisc depends on the velocity field of 
the ambient fluid, which in turn is coupled to the motion of the filament. For 
simplicity, we do not solve the full hydrodynamic problem, but instead use 
resistive force theory for slender bodies [45]. Resistive force theory amounts 
to the leading terms in an expansion in the aspect ratio a/L of slender-body 
theory, which has a nonlocal relation between force and velocity due to in- 
compressibility. To leading order, the nonlocahty can be neglected, leading to 
the local drag law, 

fvisc = -C±(V - §3 63 • V) - C||e3 63 • V, (27) 



Likewise, we take the viscous moment to be proportional to the tangential 
component of the angular velocity of the material frames, 

nivisc = -CrGs 63 • w. (28) 



The friction coefficients in Eqs. (27) and (28) are 



where rj is the viscosity of the ambient fluid, L is the total contour length of 

the rod, and a is the rod radius. For simplicity, we disregard the anisotropy 
and define ( = = (\\ = 47r?7/ log(Lo/2a), where Lq = L{t = 0). Note that 
we have kept only the leading order terms in the logarithm of the initial as- 
pect ratio. Thus, fext = — Cv- These assumptions lead to qualitative differences 
with the exact theory. For example, the neglect of hydrodynamic interactions 
implicit in the local drag approximation of Eqs. (27) and (28) will affect the 
time-dependence of the shape of the rod near self-contact points just before 
contact (see e.g. [30]). Also, the assumption of isotropy implies that the center 
of mass of a deforming closed loop remains fixed [46], whereas in the exact 
theory the center of mass can move. These limitations of the simphfied hydro- 
dynamic theory do not prevent it from capturing the essential physics of the 
phenomena we wish to study, such as the onset of buckling instabilities and 
the subsequent evolution of complex shapes. 

It is convenient for the numerical calculations of section 3.4 to write the equa- 
tions of motion in terms of position r and twist strain O3. To this end, sub- 
stitute the constitutive relation (22) into the moment balance equations (24) 
to find the force F on a cross-section 

F-^aJ5 + C(n,-.„)^x--A-, (30) 
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and the balance of the tangential components of the moment per unit length, 



C- 



ds 



(31) 



The unknown function A(s,t) occurs in Eq. (30) because the moment balance 
equation (24) does not determine the tangential component of F. Combining 
Eq. (30) with the force balance equation (23) and the expressions for fext yields 



c 



dr 
di 



so 



ds'^ ds 



dv d'^r 
^° ds ds'^ 



d_ 
ds 



A 



dv 

ds 



+ fe 



self- 



(32) 



Equation (32) has the same form as the corresponding equation for the over- 
damped dynamics of an inextensible rod, but in Eq. (32) the s-domain (length) 
depends on t. To determine A, evaluate 63 ■ dw/ds = a using Eq. (32): 



9^A d'^r d'^r . , dr d^r ^ dr 



ds'^ ds^ ds^ 



ds ds^ ds 

dr d'^r d^r 

^ ds ds"^ ds^ 



self 



(33) 



One can show that the function A in Eq. (33) is identical to the Lagrange 
multipher function of Eq. (25). 

To complete the determination of the dynamical equations, the torque balance 
(31) and compatibility relation (20) yield 



dn. 



dt 



so 



ds"^ ^ ds ds"^ dt ' 



(34) 



where the twist diffusion constant D = C/^r. Since resistive force theory 
includes the leading order terms in the expansion in a/L of the hydrodynamic 
drag force and torque, our equations are asymptotically consistent. 

The boundary conditions which accompany Eqs. (32-34) depend on the situa- 
tion. For a filament with free ends, the appropriate conditions are M = and 
F = 0. For a closed loop, the variables r, 17, w, F, and M must be periodic 
in So with period 2t:Rq. 



2. 5 Change of Basis 



We close this section by returning to the claim of section 2.2 that the spon- 
taneous twist To does not affect the shape of an inextensible rod with circular 
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cross-section and vanishing spontaneous curvature kq = 0. Consider such a rod 
with n = ToGs in the stress-free state. We can ehminate Tq from the problem 
using the freedom to redefine the material frame. If 



e[ = Gi cos — §2 sin i 
§2 = §1 sin 4> + e2 cos i 



(35) 
(36) 



then 



Q'l — Qi cos (f) — Q2 sin 
^2 = —^1 sin -|- Jl2 cos (f) 



ds 



and 



(37) 
(38) 

(39) 



cu'i — cui cos (f) — UJ2 sin 
UJ2 — —0J\ sin -|- a;2 cos 



so 



Under this change of basis, fi'j^e'^ -|- ^22^2 = ^iGi + ^2^2 for any 
= — STo fixes = in the stress-free state. Thus, 



(40) 
(41) 

(42) 

!). Choosing 
(43) 



the parameter Tq has been eliminated from the constitutive relation. Note that 
our argument up to this point holds for a growing rod as well. 

Now consider the effect of the transformation (35-36) on the compatibihty 
relations. Once again, even in the presence of exponential growth, the com- 
patibility equations take the same form, e.g. 

= — ^ - + Vl[uj'2 - 9!^uj[. (44) 




For an inextensible rod, uj'^ — uj^, since s — Sq ii a — Therefore, the 
change of basis (35-36) does not affect the rotational drag or translational 

drag equations, and we conclude that Tq is not a physical parameter for an 
inextensible rod with circular cross-section and no spontaneous curvature. 
However, we expect the opposite conclusion for a growing rod, since we saw 
in section 2.2 that Tq has physical meaning. In fact, once Tq is eliminated from 
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the constitutive relation using Eqs. (35-36), a new To-dependent term appears 
in the torque balance equation (31): 



C 



(45) 



ds 



Therefore, Tq cannot be eliminated from the equations of motion for this mode 
of growth. 



3 The Growing Elastic Loop 

In this section wc treat the problem of a growing elastic loop with preferred 
twist Tq. As described in the introduction, this example is motivated by the 
growth of B. subtilis filaments from a spore. Sometimes the ends of the grow- 
ing filament stick to the spore coat, leading to a closed loop [18]. The filament 
lengthens, depleting the twist, but the closed geometry prevents the rotation 
of cross-sections normally seen in unconstrained filaments. Thus, twist stress 
builds up, and the filament eventually writhes and coils. Although the model 
does not address the writhing and coiling of unconstrained filaments, it dis- 
plays some of the features exhibited by B. subtilis loops. 



3.1 Buckling and Writhing Instabilities 

To study the stability of an exponentially growing circular loop with preferred 
twist To, we begin with the unperturbed solution. The unperturbed loop lies 
in the z — plane and has radius R — Roexp{at). Since each material point 
on the filament moves radially outward with fixed z in unperturbed growth, 
the angular velocity vanishes, lj^^^ = 0, where we use the superscript, (0), to 
denote the unperturbed value. The bending part of the energy (25) decreases 
exponentially as the loop grows because k^^^ — 1/R. However, the twist en- 
ergy density increases because the closed geometry prevents the cross-sections 
from rotating with the rate asro required to attain the twist state of zero 
energy. If we assume for simplicity that Qf\s,t = 0) = tq, then the expo- 
nentially increasing length leads to an exponentially decreasing twist density, 
ilf\s,t) = Toexp(— at). To summarize, the moment on a cross-section takes 
the form 




(46) 
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Note that it is much more convenient to express M in terms of the cyhn- 
drical coordinate unit vectors {z, p, (p} instead of the material frame vectors 
{e^'''*, 62°'', ef^}, since e^^"* and e'^^ continuously rotate about eg^'' as s increases. 
A choice of frame which does not rotate about the tangent vector as arclength 
increases is a natural frame [47] . We will use the natural frame extensively in 
the linear stability analysis of section 3.3. 

Since the filament is simultaneously bent and twisted, moment balance (24) 
implies a component of force in the 2;-direction, whereas force balance (23) 
leads to a tangential component proportional to the growth rate a: 

= -C(7i?V + — (e-* -l)z. (47) 



The tangential force on the cross-section is compressive and grows exponen- 
tially in time, eventually leading to an Euler-like buckling instability when 
CaR^ « A/R^, or R = Ri^ [^/(Ca)]^/^ (c/. [30,29,31]). The Euler buckling 
time scales as ~ log[A/ {(aR^)] (when A/{(aRQ) < 1, the loop begins 
to buckle at i = 0). In section 3.3, we will refine this estimate using our linear 
stability calculation and see that the appropriate buckling time at small a 
actually scales as cr~^, since in this regime, a growing perturbation does not 
become noticeable until long after the perturbation begins to grow. However, 
our numerical calculations of section 3.4 reveal that the correct picture is even 
more complicated: for sufficiently small TqRq, a secondary instability arises 
and grows before the linear instability has significant amplitude. In addition 
to the Euler buckling, there is also a writhing instability. The magnitude of 
the twist moment increases as |Cro(l — exp{—at))\ as the filament lengthens, 
leading to a writhing instability (like that of a twisted ring [48,49,50]) when 
Cto{1 — exp(— crt)) pa A/{27rR), or % ~ cr~Mog[l + l/(ro-Ro)], (assuming 
A Ki C). Note that the critical time for writhing is % ~ ^og[l / (RoTq)] 
for RqTq <S 1, and % ~ 1/(c-RoTo) for RqTq ^ 1- The sense of rotation of 
the cross-sections of unconstrained filaments determines the handedness of the 
coils that form after the instability of the ring: positive tq (counter-clockwise 
rotation when viewed from the direction of increasing s) leads to right-handed 
plectonemic braids. Note that drag is the ultimate cause for the Euler buck- 
ling instability; once growth ceases, the buckled filament relaxes back to the 
unperturbed circular shape. Since the writhing instability arises not from drag 
but instead from the frustration of growth-induced twist stress, the braided 
post-instability shape remains after growth ceases. For small enough TqRq, we 
will see in section 3.4 that the small intrinsic twist biases the Euler buckhng, 
leading to solenoidal shapes, which relax to plectonemes if growth halts af- 
ter a sufficiently long time. Observations of the growing fibers suggest that 
writhing is the dominant mechanism in the instability of a closed loop [18]. 
Using A = 10-^2 erg-cm [28], a = 2x 10"^ s-^ L = lO'^ cm, a = 3 x 10"^ 
cm, and ( = fa IQ-^ erg-s/cm^ (see Eqn. (29)) leads to [A/(C(t)]^/^ 150 
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//m, which is much larger than the observed critical radius and implies that 
growth-induced Euler buckling is not important. Thus, we can use the criti- 
cal radius of the twisted ring to estimate Tq; assuming tq is comparable to or 
smaller than I/-R0, where Rq is the initial radius, and using i?2 ~ 2 x 10~^ 
cm as the experimentally observed critical radius [18], we find tq ~ 10^ cm~^. 
This twist rate corresponds to a few turns per cell. It is intriguing to note that 
the corresponding length scale is close to the pitch of helical filaments of mbl, 
a recently discovered protein which resides near the cell wall and plays a role 
in maintaining the shape of B. subtilis cells [51]. 



3.2 The Natural Frame 



Before studying the evolution of small perturbations of the growing circular 
shape, we consider the choice of representation. Since only the shape and twist 
are of interest, it is convenient to use an intrinsic representation, such as the 
material frame. However, we saw in the solution of the unperturbed problem 
of section 3.1 that the natural frame leads to further simplification. The ad- 
vantages of the natural frame over the material frame are even greater for 
the linear stability analysis, and more generally for the full nonlinear problem 
[42,52]. 

In a natural orthonormal frame {fii, n2, 63} the instantaneous rate of rotation 
of 111 about 63 is zero [42,47,52]: 

...f^O. (48, 



A rotation of {ni,n2} about 63 by uniform (arclength-independent) angle 
leaves the condition (48) invariant; every space curve has a family of nat- 
ural frames, the members of which are related to each other by rotation 
through a uniform angle. To construct a natural frame from the material 
frame {61,62,63}, rotate the material frame at s by minus the accumulated 
rotation angle ■»? = Jq ds'Qs: 



hi — cos i^ei — sin i?e2 (49) 
n2 = sin i?ei -|- cos ■J?e2 . ( 50 ) 

The natural frame is nonlocal in the sense that deformations of the filament 
center line in the region s' < s affect the natural frame at s. 

Our formulas can be further simplified with complex notation. For example. 
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if e = ni + in2, then Eqs. (49-50) become 

e = (ei + ie2) e^'' . (51) 



The vector e is a complex normal vector. Note that e • e = 0; therefore, 
e • de/ds — 0. Furthermore, the defining property fii • 8x12/ ds — implies 
e* • de/ds — 0. Thus, de/ds is proportional to 63. In analogy with the Prenet- 
Serret equation dh/ds — —Kt, we define the complex curvature ^ via 

I = -*e3, (62) 



where 



= (-iOi + O2) e^'^ . (54) 

Also, the rate of change along s of the unit tangent vector is the complex 
curvature times the complex normal vectors, 

f =l(*e- + *-e). (65) 



Note that rotation of a natural frame about 63 by a uniform angle leads to a 

constant shift in the phase of ^. For example, the natural frame arising from 
the construction (49-50) applied to the Frenet-Serret normal and binormal 
has 

e = (n + ib) exp i / rds , (56) 



since torsion is the rate at which fi and b twist around the tangent vector. 
The corresponding complex curvature is 

*FS = «;exp (i /rds); (57) 



the ratio of ^ and ^^ps is a constant phase. 

To complete the specification of the kinematics of the natural frame, we in- 
troduce the complex angular velocity describing the rate of change with time 
of the unit tangent vector: 
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n = e- 



de3 

dt 

- + a;2 I e 



(58) 
(59) 



In terms of the natural frame variables, the compatibihty relations (18-20) 
become 



9^ 


dn 


'dt 


ds 

so 




_ dujs 


dt 


ds 



s 

(7* - i^a;3(0) + / ds'lm(^**n 



- aQa + Imf 



so 



(60) 
(61) 



The integral in Eq. (60) reflects the nonlocality of the natural frame, and arises 
from the temporal rate of change of 



dd_ 
'dt 



so 



d_ 
di 

d_ 
di 



soe" 



-I 



n3{s',t) ds' 

so 

so 

J 1^3(4 e'^*,t)e'^* ds'o 
ds'; 



so 

dQ 



dt 



+ (7^3 



(62) 

(63) 
(64) 



Eq. (60) follows from Eqs. (64), (18-20), and the identity fliu!2 — f^2'^i = 
lm(^*n). 



In the natural frame variables, the moment M obeys 



(65) 



Just as in section 2.4, moment balance (24) determines the perpendicular 
component of the force on a cross-section 



F = F||e3 + ^(e*F^ + eF:); 



(66) 
(67) 



Recall that Fy is not determined by moment balance since only the perpen- 
dicular components of F enter Eq. (24); F\\ is determined by the condition 
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63 • dv/ds — a. 



Defining the tangential and perpendicular components f\\ and f± of the force 
per unit length 



dF 

ds 



63/11 + ^(e7± + e/l), 



(68) 



it follows from Eqs. (52), (55), (66), and (67) that 



(69) 



and 



52^ d 



^ - To 



(70) 



Note /|| = d{F\\ +A\^f\'^/2)/ds = -dA/ds; in terms of /y and f± the condition 
■ dv / ds — a becomes 



^ ds 2-' 2-'^ 



(71) 



Finally, the force per unit length determines the complex angular velocity 
through the relations = dF/ds and e • dv /ds = U, or 



df± 



(72) 



Likewise, the tangential component of the angular velocity is given by Eq. (31). 
In summary, the equations of motion for the filament in the natural frame are 
Eqs. (60-61), (69-72), and (31). 



3.3 Linear Stability Analysis 



We now return to the stability analysis of a growing circular ring. We write 
the shape as r = (i? + r'^^^)p + z'^^^z and work to first order in r^^^ and z'^^^ 
It is convenient to find the equations of motion for ^^^^^ and ^3^^ first, and 
then express these equations in terms of r^^^ and z^^\ To this end, we choose 
{hf\ np} = {-p, z} to make = 1/R. Expanding (57) to first order leads 
to 

^ = «;(0)+«;(l)+i«;(0) Jr^^Ms, (73) 
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since the unperturbed loop has zero torsion, r*^°) = 0. To express k and r in 
terms of r^^^ and 2;^^^ use the Frenet-Serret equations (8-9) and s — R(j) (to 
leading order) to find 

+ (74) 



and 



Thus, * = + + i77, with 

e = -;^(r(^)+rS) (76) 

^=:^(^^'^+4'i)- (77) 

As described in section 3.1, the unperturbed angular velocities vanish: 11*^°) = 
0, cl;3'^'' = 0; furthermore, Eq. (47) together with (69-70) imply = F'i^^^ / R = 

—(aR and fu^^ = 0. Expanding the equations of motion to first order, we find 



dt 



and) 



ds 

so 



s 

_ ^^(1) _ i*(o)a;(^)(0) + i^(°) J d/Im(*(°)*n(i)) (78) 



dt 



D^0L _ an^) + Im(^(°)*n(^)) (79) 



so 



= ^_Re(/f -Re(/l^)*(o)*) (80) 

^nw = ^ + *(o)/W, (81) 

OS " 

where fj^^ and /j^"* are determined by expanding (69-70) to first order. Inspec- 
tion of (78-81) and (69-70) reveals that r^^^ = r^^^ cos(n0), z^^^ = z^^^ sin(n0), 
Fy"^^ = F|j cos(n0), and = cos{n(j)) , with n a positive integer. This 

choice of origin for s and Eq. (31) imply 0)3^^ (0) — 0. The perturbation r^) cos 
corresponds to a translation of the ring in the z — plane; likewise, the per- 
turbation sin (j) corresponds to a rotation of the ring about an axis in the 
z = plane. Thus, the n = 1 perturbations are rigid motions, leading to no 
change in curvature: C, = i] = (see Eqs. (76-77j). Inserting the perturba- 
tions into (78,79,81) and using (80) to ehminate yields a linear system of 
differential equations for r^^\ 5^^, and Q3: 



(1) 
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5(1) 



n2 + 1 + 1 



£(1) 



c 

c 



(82) 



(0) 
3 



To)n(n2 - l)fW + 



3t ■ 



C 



+ 



an 



n^n^ - 1) 



(83) 



— a 



nil — n 



(84) 



Equations (82-83) hold for n 7^ 1, and equation (84) holds for all integer 
n > 1. For n — 1, the shape drops out of Eq. (84) since the n — 1 shape 
perturbations are rigid motions. Note that Tq enters (82-84) in the combination 
^if^ — Tq = ro(exp(— (xt) — 1); tq drops out of the equations when cr = 0, in 
accord with the general arguments of section 2.5. 

We can simphfy the hnear system (82-84) by exploiting the large ratio of 
relaxation time scales for twisting and bending modes: 



''bend 



> t 



twist ) 



(85) 



for R ^ 10 /im and the B. subtilus parameters of section 3.1, tbend = C-R^/^ ~ 
10~^ s and itwist = R^/D — (^,R'^/C 10~^ s. As time passes and the filament 
lengthens, tbcnd and ttwist increase exponentially, but tbcnd ^ ^twist for all 
time. Therefore, twist perturbations Q relax immediately, and (82-84) may 
be simplified by setting Cl to zero. Using the initial radius for the bending 
relaxation time (tbend = C-^oM) assuming C/A — 1 ior simplicity, (82- 
83) reduce to 

q = ^q, (86) 



where 



q - ( ^(1) ) (87) 



and 



-"1+3 



•TQ-Roe 



^bend 



-(e-'^*-l)n(n2 - 1) 



^bend 

or? — 



^bend 



at 



n^iri^ 



1) 
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Since L depends on time, the system (86) is non-autonomous, and classical 
modal analysis [53] docs not apply. Note also that L is not a normal operator: 
[L, L^^] 7^ 0. Therefore, the eigenvectors of L are not perpendicular, which in 
general signals the possibility of transient or algebraic growth of perturba- 
tions [54]. However, L is only weakly non-normal, since the eigenvectors are 
almost perpendicular. A similar conclusion applies to the problem of a ring 
with twist but no growth {a = 0) in a viscous fluid (c/. [50,43]). Therefore, 
we do not expect the phenomenon of transient growth of perturbations in 
our problem. Nevertheless, the methods developed to study non-normal linear 
problems are well-suited to non-autonomous problems such as (86) [55] . 

We will characterize the growth of perturbations by the amplification of the 
magnitude of q(0). The optimal amplification G{t) is defined by maximizing 
this factor over all initial conditions: 

G(t)^ max 44^4^. (89) 
^ ^ q(o) q(0) • q(0) ^ ^ 



(A more physical choice for the optimal growth rate would be to use the 
bending energy to second order in f^^^ and z^^^ instead of q • q; it turns out 
that either choice yields essentially the same G{t).) To compute G{t), recast 
Eqn. (86) as an equation for the propagator matrix B: 

B = LB, (90) 



where 6(0) = / and q(t) = 6(t)q(0). Thus, the optimal growth factor is the 
Rayleigh quotient 

^ ' im q(0) ■ q(0) ^ ' 



or G(t) = A+. where A_|_ is the largest eigenvalue of B'^{t)B(t) [56]. When 
/- is a t- independent normal matrix, then A_|_ = exp(2A_|_t), where A_|_ is the 
largest eigenvalue of L. Note that since G{t) is maximized at each t over all 
initial conditions q(0), the maximum amplitudes at two different times may 
correspond to two different initial conditions q(0). We computed G{t) by using 
standard Runge-Kutta techniques to solve for B{t) [57]. 

Inspection of the diagonal components of L (88) reveals that for sufficiently 
rapid growth, cr > (Tc = n^(n^ — l)^/[(n^-|-3)ibend], the loop deforms away from 
its circular shape as soon as it begins to grow. When the rate of growth of the 
ring is sufficiently slow, a < a^, bending stiffness stabilizes the circular shape 
for crt -C 1 and perturbations decay roughly as exp(— n^t/tbend)- Thus, the 
growth factor G decays extremely rapidly with increasing n at early times. As 



23 




JQ I I I I I I 

0.5 1 1.5 2 2.5 3 

ct 



Fig. 5. Optimal growth factors for atbend = 0.5 (< (Tcibend)- Solid line: -Ro'^o = 0. 
Dashed line: RqTo = 5. Dotted line: i?oTo = 10. Note that the intermediate regime 
of rapid growth of G{t) is most apparent for RoTq = 10, especially for n = 3. 

time passes, the loop lengthens and eventually buckles, with the nature of the 
buckling dependent on the magnitude of tq. For tqRq <C 1, the distortion is 
the three-dimensional analog of the in-planc Eulcr buckling studied by Shelley 
and Ueda [29,30] and Drasdo [31]. For tqRq ^ 1, the off-diagonal elements of 
L are large [see (88)], writhing dominates the nature of the initial distortion, 
and G{t) increases roughly as exp(2ro-Ro'^^V^bend) in the intermediate regime 
at <Ri 1. In the late-time regime at » 1, G{t) oc exp(2n^(Ti) for any value of 
toRq or t]j. These results are summarized in Fig. 5. 

There are two different times which may be chosen to represent the time of 
the instability. The estimate for the buckhng time of section 3.1 amounts to 
the time ti at which G{t) reaches its minimum value and starts to increase. 
However, since the rate of decay of G{t) in the stable period can be very 
different from the rate of growth of G{t) in the unstable period, the time at 
which a perturbation becomes noticeable may be significantly greater than ti. 
Thus, it is natural to define the time for the onset of the instabihty to be the 
time t2 at which the amplitude of the optimal perturbation regains its initial 
value: G{t2) = 1. The graphs of Fig. 5 suggest that t2 and ti are comparable 
whenever -RqTo is large enough for the intermediate growth regime of rapid 
growth discussed above to be present. However, if RqTq is small enough that 
this intermediate regime is absent, then t2 will be much greater than ti when 
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Fig. 6. Time t2 for optimal perturbation to regain initial amplitude as a function 
of ii!oTo for atbend = 0.001 (solid line), attend = 0.01 (dashed line), aibend = 0.1 
(dotted line), crtbend = 0.5 (dashed-dotted line), and crtbend = 1 (solid line with 
small dots). The large filled dots arc the results of numerical simulations of the full 
nonlinear equations with aibend = 0.01 and the open squares are with utbend = 0.1; 
see section 3.4. 

the extensional growth rate is slow, (xtbcnd ^ 1- Figure 6 shows how dramatic 
this difference can be. For RqTo -C 1 and unbend 1, 



whereas ti oc 1/a. When RqTq ^ 1, both ti and t2 scale as l/o" (see Fig. 6). 
Thus, for small (xtbend, there is a sharp transition in the onset time t2 as a 
function of RqTq. In section 3.4 we will see how this prediction of the linear 
theory captures the early-time dynamics for RqTo ^ 1, but that nonlinearities 
intervene before t — t2 for RqTq <S 1. 

The curves for the onset time t2 of Fig. 6 were computed from the linearized 
equations (86) using the adiabatic theorem. If at bend ^ 1, and if L were 
normal, then the adiabatic theorem [58] would imply that 




(92) 



aH^end 4(n4 + 3) ' 




(93) 
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Fig. 7. Optimal growth factors for tqRq = 10 and crtbend = 1-0, 0.5, 0.1. 

where v±v± are the dyads formed from the eigenvectors v±(t) of L(t). Since 
L is not normal, Eqn. (93) is in error by an amount governed by v+ • v_, 
which is never more than about 0.1 and is often much smaller. Note also that 
if crtbend ^ 1) then the off-diagonal elements of L are small compared to the 
diagonal elements, causing the equations for f^^^ and z^^^ to decouple and 
leading to 

G{t) ^ exp J A+(t')dt'j , (94) 



where A_|_(t) is the largest eigenvalue of L{t). Thus, in both extremes crtbend ^ 1 
and attend ^ 1; G{t) ^ exp(2 / A+dt'). This result is especially useful in the 
limit of small growth rate utbend <S 1, since the rapid relaxation and growth 
of bending modes makes it difficult to solve for B numerically. 

Figure 7 shows the optimal growth factors for the first two modes (n = 2 and 
n = 3) for crtbend — 1-0, 0.5, 0.1 and RoTq — 10. For both the intermediate 
writhing regime and the large-crt asymptotic regime, the rate of increase of the 
growth factor increases with mode number n. As utbcnd increases, the time at 
which the growth factor for the n = 3 mode overtakes that of the n = 2 mode 
decreases because the instability of each mode occurs at earlier times. For the 
larger values of crtbend) the n — 3 mode overtakes the n — 2 mode before the 
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Fig. 8. Shape of the growing loop for tqRq = 10.0, 0,3{t = 0) = 8.3, and atbend = 0.1 
at (1) t/Wd = 0.03, (2) t/tbend = 2.41, (3) t/Wd = 3.22, and (4) t/Wd = 3.37. 

amplitudes have grown large enough for nonlinearities to come into play. Thus, 
we expect to see double-stranded plectonemic braids with two hairpin turns 
for small crtbend, and braids with three or more hairpin turns when crtbend is 
large. The numerical computations of section 3.4 confirm these expectations. 

3.4 Numerical Solution of the Nonlinear Equations 

We solved the closed set of equations (32,33,34) using a pseudospectral method 
[59] for the backbone dynamics (32), direct integration of (33) at each time 
step using finite differences to find A, and a Crank-Nicholson routine for the 
twist dynamics (34). We used initial conditions such that the backbone of the 
loop is perturbed from circular shape with i?o = 1 by a few small-amplitude 
modes (n = 2-5). Depending on the simulation, ^3 ranged from somewhat less 
than To to tq. Fig. 8 shows a time series of the shape of the growing loop with 
tqRq = 10, Q^it = 0) = To, and atbend = 0.1. For early times (Fig. 8-1), the 
circular loop is stable and perturbations decay. As the loop grows, R increases 
exponentially and fl^ decreases. At a critical value of fl^ and R (Fig. 8-2), the 
loop begins to buckle and wrap about itself. For sufficiently large tqRq, the 
loop takes on the conformation of a plectoneme, initially forming a figure-eight 
structure (Fig. 8-3) and then wrapping into a braided form (Fig. 8-4). Fig. 9 
shows the twist energy, / C^Q-^ — tq)^/2 ds, and bend energy, / ds, for 

the growing loop depicted in Fig. 8. Note that the total energy is not fixed 
in our model since growth acts to inject energy into the system. At point 
1 in Fig. 9, growth along the filament axis leads to a decrease of twist in 
time and thus an increase in the twist energy. At the same time, backbone 
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Fig. 9. Semi-log plot of the bend (solid line) and twist (dashed line) energy vs. time. 
Numbers correspond to shapes in Fig. 8 and label important times during growth 
(See text for further explanation). 

perturbations die away and the curvature decreases exponentially, leading to a 
decrease in the bend energy. At point 2, the circular loop becomes observably 
unstable. The bend energy begins increasing as perturbations in the filament 
grow, and the twist energy continues increasing (See Fig. 9). At the inflection 
point, point 3, the filament forms a figure-eight pattern. Note that a figure- 
eight shape of a closed loop which is not growing is a minimum of the total 
energy for a range of twist. At later times (such as point 4), the filament wraps 
into a braided structure. The bend energy increases as more braids are added. 
The twist energy also increases; however, writhing motions act to decrease 
the twisting stress imposed by growth, leading to a twist energy that grows 
sub-exponentially (see Fig. 9). 

The dynamic equations were solved for a range of tqRq and crtbend- When 
To-Ro was larger than about 4-10 (the actual value depended on the value 
of (xtbcnd), "we found plectoncmcs. Smaller values of toRq yielded solenoids. 
Figure 10 shows how the morphology of the loop depends on the value of TqRq 
at fixed attend- In the region with plectonemes, the shapes form as a result 
of a writhing instability; the Euler-like buckling instability [29,30,31] plays 
little role. In particular, the braided shapes remain after growth ceases. Thus, 
the shapes are qualitatively similar to the minimizcrs of the elastic energy 
without growth. However, as discussed in section 3.3, the rate of growth affects 
the shape since the number of branches increases with attend- In the region 
with solenoids, the Euler-like buckling instability comes into play since the 
small value of tq-Ro delays the onset of the writhing instability. Solenoids are 
not minimizers of the elastic energy without growth; when growth ceases the 
solenoids relax. Therefore, the solenoids are the three-dimensional analogs of 
the two-dimensional shapes of references [29,30,31]. The plectonemic, three- 
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Fig. 10. Morphology of the growing loop at different values of tqRq. Solenoidal 
morphologies occur when roi?o < 2. Branched plectonemes are observed for 
2 < tq-Ro ^ 3 and unbranched plectonemes are observed when roi?o > 4. 
c^bend = 0.01 for all figures. 

armed, and solenoidal morphologies we obtain are similar to many of the 
supercoiled patterns that are observed in B. subtilis; in the conclusion we 
discuss the relation between our results and the experiments. 



In these calculations, we defined the time to the onset of the instability as the 
point where the bend energy begins increasing (point 2 in Fig. 9). This time 
should correspond roughly to t2 defined in section 3.3, since the perturbations 
must be noticeably large before they affect the bending energy. Just as in the 
linear stability analysis, increasing tq-Ro or a reduces the time to onset of the 
instability. For values of tq-Ro > 10, the time to onset found numerically was 
in quantitative agreement with the linear stability analysis (see Fig. 6). At 
smaller values of tq-Ro the onset time from the simulations was earlier than 
that predicted by the linear analysis, with the deviation getting larger for 
decreasing atbend- The deviation is due to a secondary instability: before the 
growing n = 2 mode becomes observably large, it is overtaken by a higher 
order mode which quickly dominates the shape of the loop. The value of the 
mode that dominates depends on the value of crtbend and tqRq. As tqRq is 
decreased at constant crtbend; the mode that dominates increases. As crtbend 
is increased at constant Tq-Ro) the mode that dominates decreases. Figure 11 
shows the phase diagram and examples of plectonemes and solenoids. As a 
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Fig. 11. Qualitative phase diagram depicting the dominant mode behavior as a 
function of crtbend and tqRq. For tqRq greater than roughly 10, the linear stability 
analysis predictions are valid and the n = 2 and n = 3 modes dominate the pattern 
formation. At lower values of tqRq, a secondary instability drives pattern selection 
and higher order modes dominates. Insets show the morphology of the loop in 
different regions of the diagram. 

check, we verified the secondary instability by discretizing the full nonlinear 
equations (32,33,34), and integrating with the MATLAB routine odel5s (a 
variable order, variable time step method for stiff problems). 



4 Conclusions 



The three essential elements of any growth process appear in our model for 
B. subtilis fibers in a simple form: the mode of growth is exponential exten- 
sion with rotation, the material properties are described by the bending and 
twisting elasticity of a slender filament, and the interaction with the external 
environment is governed by resistive-force theory. Thus, our model is ideal for 
illustrating the basic phenomena of the physics of growth. We have assumed 
that the bacterial filaments can be treated as perfectly elastic and that the 
growth rates are uniform and independent of stress. The degree to which these 
simplifying assumptions hold remains an outstanding experimental problem. 
Furthermore, our model does not lead to plectonemic structures for the case 
of bacterial fibers with free ends for the measured values of the growth param- 
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eters and elastic moduli. Nevertheless, we can draw some general conclusions 
from our model and use our results to suggest new mechanisms for pattern 
formation in the presence of biological growth. 

First, even though the ultimate microscopic mechanism for the supercoiled 
patterns of B. subtilis fibers is unknown, the blocked rotation mechanism we 
study here must play a role in the supercoiling of closed fibers. For exam- 
ple, since the closed loops form supercoils at lengths which are much shorter 
than the length at which the open fibers supercoil [19], we expect that the 
blocked rotation mechanism dominates over whatever mechanism causes the 
open fibers to supercoil. Thus, we can directly compare our calculations with 
the experiments on the supercoiling of closed filaments. Recalling the param- 
eters of section 3.1, wc estimate a%end 1, so that we expect unbranched 
plectonemes when TqRq > 4, branched plectonemes when 2 < tq-Ro ^ 4, 
and solenoids when tq-Ro ^ 2. Typical observations of the filaments yield 
unbranched plectonemes, but branched plectonemes and solenoids also arise, 
depending on the growth medium. Using the observed buckling radius for an 
unbranched plectoneme, in section 3.1 we were able to estimate Tq ~ 10^ cm~^, 
suggesting the presence of a structure in the cell wall with pitch P ^ 10~^ cm. 
The observations of helical actin-like polymers in the cell wall with compara- 
ble pitch [51] support our estimate, and further suggests a starting point for 
a theory for the microscopic mechanism of the supercoiling. 

The second major conclusion of our work is the dynamical nature of the pat- 
tern selection. For large TqRo, the plectonemes remain once growth ceases, and 
are qualitatively similar to the minimizcrs of the elastic energy, although the 
rate of growth plays an important role in determining the shape. For small 
To-Roi the solenoids arc transient structures which relax away when growth 
halts. Since the curvature of the bacterial fibers can become permanent [26], 
the solenoidal shapes may act as a template for patterns which remain in the 
absence of growth. This mechanism of pattern formation via the hardening 
of transient structures formed from the interplay of flexibility and external 
friction may apply to other biological systems, such as those studied in [31]. 

To help justify, refine, or rule out our model, we suggest three basic experi- 
ments. First, the shape of a growing loop as function of time should be mea- 
sured precisely enough to compare with our theory. Although our numerical 
results (Fig. 8) arc qualitatively similar to the experimentally observed shapes 
(Fig. 1), the lack of detailed information about the evolution in time of a sin- 
gle loop in [18] prevents a stringent test of our theory. Second, the change in 
elastic properties with time during growth should be quantitatively measured. 
For example, to what degree and how long must a fiber be bent to develop a 
permanent curvature? Finally, the nature of the twist stress in a growing fiber 
with free ends should be determined. Is there a twist moment on the cross 
sections of the growing fibers with free ends, leading to a writhing mechanism 
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qualitatively similar to the blocked rotation mechanism for closed loops, or is 
the mechanism for open fibers completely different? Future progress toward 
understanding the pattern formation of B. subtilis, both at the microscopic 
level of the structure of the cell wall and the more macroscopic level treated 
here, depends critically on new experiments such as these. 
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